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ABSTRACT 

We present the first simulations of evolving, strongly twisted magnetar magnetospheres. Slow shear- 
ing of the magnetar crust is seen to lead to a series of magnetospheric expansion and reconnection 
events, corresponding to X-ray flares and bursts. The axisymmetric simulations include rotation of 
the neutron star and the magnetic wind through the light cylinder. We study how the increasing twist 
affects the spindown rate of the star, finding that a dramatic increase in spindown occurs. Particu- 
larly spectacular are explosive events caused by the sudden opening of large amounts of overtwisted 
magnetic flux, which may be associated with the observed giant flares. These events are accompanied 
by a short period of ultra-strong spindown, resulting in an abrupt increase in spin period, such as was 
observed in the giant flare of SGR 1900+14. 

Subject headings: magnetic fields — plasmas — relativistic processes — pulsars: general — 
magnetohydrodynamics 



1. INTRODUCTION 

Magnetars are neutron stars with luminosity 
powered by the di s sipati on of magnetic energy 
(|Duncan k. Thompson! I1992D . Their magnetospheres 
evolve on timescales much shorter than the age of the 
star. The magnetosphere is attached to the stellar crust 
and is thought to be deformed by crustal motions (either 
gradual or catastrophic). 

The crust is practically incompressible, and its mo- 
tion must be pure shear. Shearing of the magneto- 
spheric footpoints generates external twists (V x B / 0) 
and hence excites external electric currents. Ohmic 
dissipa tion of these currents generates non-thermal ac - 
tivity ([Thompson et all 120001 120021 : lBeloborodovll2012h . 
A quasi-steady twist implanted in the magnetosphere 
is erased on the Ohmic timescale ~ 1 yr, as it 
gradually dissipates via a contin u al electric di s charg e 
(Beloborodov fc Thompsonl [20071 : iBeloborodovl 120090 . 
Thus, the crustal shearing creating the twist must be 
sufficiently fast, with vorticity ui c <; 1 rad yr . 

Observed magnetars rotate with periods P = 2 — 12 s 
and are gradually spun down, on a timescal e compara- 
ble t o the age of the star, t ~ 10 4 yr (e.g. Me reghettH 
2008). The spindown torque acting on the star is con- 
trolled by the open magnetic flux passing through the 
light cylinder of radius i?LC = cP/2ir. As the magneto- 
sphere is twist ed, its energy density grows and it tends 
to inflate (e.g. lWolfsonlll995h : this opens more magnetic 
flux, increasing the spindown torque. Thus, the magnetic 
twists, besides generating non-thermal emission, can also 
account for the observed temporal variations in spindown 
rate P. 

For example, an outburst from XTE J181 0-197 was 
followed by a chan ge in P of more than a factor of three 
over nine months (Camilo et al.1 12007T ). The spindown 
rates of SGRs 1900+14 and 1806-20 have been observed 
to vary by a facto r of four over timescales of months 
(jWoods et aU2 002). The 27 August 1998 giant flare from 



SGR 1900+14 was coincident with a fractional increas e 
in the spin period of AP/P = 10" 4 (jWoods et al.|[l999l) . 
Because of an 80 day gap in observations before the flare, 
the behavior of P is unknown. It is unclear if there was 
a gradual change in spindown over this period, if AP 
resulted from a brief and dramatic increase in magnetic 
torque, or if it was a result of a sudd en change inside 
the n eutron star, a sort of "anti-glitch" ([Thompson et all 
[2003). 

In a simple axisymmetric model, the twist amplitude 
tp is measured by the azimuthal angle between the mag- 
netospheric footpoints. Twists of small amplitude )/i < 1 
change P by a small fraction ~ ip 2 (|Beloborodovi r2009) . 
Significant changes in P are associated with large ip ~ 1. 
On the other hand, when ip exceeds a critical value 
~ 1, the twisted magnetic equilibrium loses stability (e.g. 
lUzden skv 2002). The instability is similar to coronal 
mass ejections from the solar corona — part of the mag- 
netic energy is ejected, and the twist amplitude is re- 
duced. 

The purpose of this paper is to study the response of 
the magnetosphere of a rotating star to strong twisting. 
We present the first numerical simulations of this prob- 
lem, in its simplest, axially symmetric version. We con- 
sider the aligned rotator (a dipole magnetosphere whose 
axis is parallel to the spin axis of the star); the twisting 
motion of the stellar crust is modeled as a slow rota- 
tion of either one polar cap, or a band of latitudes, with 
respect to the rest of the star. 

The magnetospheric field lines are conveniently labeled 
with the poloidal flux function /; for the twisted dipole 
configuration / = 27r/isin 2 0/r+, where fj, is the magnetic 
dipole moment of the star, r* is the stellar radius, and 
8 is the colatitude of a field line's northern footpoint. 
The gradual pump ing of the twist is described by (see 
IBeloborodovl[20Tl 

W)=Wc(/) + 27rc|| =«(/). (1) 
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Here w c = <fi n — 4>s is the differential angular velocity of 
the moving crust ( "n" and "s" refer to the northern and 
southern footpoints of the field line), and $ is the voltage 
between the footpoints. The resistive term 2wcd^/df is 
negligible if ui c 3> lrad yr" 1 . When this term is impor- 
tant, it can temporarily amplify th e twist near the axi s 
while reducing it everywhere else (|Beloborodovl 120091 ). 
The approximate effect of the resistive term can be incor- 
porated into an ideal calculation by attributing the entire 
effective uj to footpoint shearing. This is the approach we 
take in our simulations, where the magnetosphere is de- 
scribed as ideal and force- free (except in current sheets). 

Our simulations are novel in at least two re- 
spects. Firstly, previous time-dependent numerical in- 
vestigations of sh eared magnetospheres have used non- 
relativistic MHD (|Mikic fc Linkerlll994t lAntiochos et al.l 
1999); we solve the equations of force- free electrodynam- 
ics, the appropriate limit for plasmas in ultra-strong 
magnetic fields. In this regime, the plasma can be 
approximated as massless. Secondly, our simulations 
include stellar rotation, with the light cylinder radius 
within the computational domain. 

2. NUMERICAL SIMULATIONS 

Our simulations are performed with PHAEDRA, a 
para llel pseudospectral c ode for force-free electrodynam- 
ics (jParfrev et al.l I2012D . The equations solved are 
Maxwell's equations with a non-linear current density J, 
which incorporates the plasma response: 



d t E 



J = 



d t B = -V x E, 
B-VxB-E-\7xE 



V x B - 4tt J, 
E x B 



AttB 2 
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(2) 
(3) 



These equations ensure that the Lorentz force vanishes, 
(V • E)E/4:'K + J x B — 0, and that E is perpendic- 
ular to both B and J, implying zero Ohmic dissipa- 
tion. Motion of the stellar surface is effected by impos- 
ing an electric field at the inner boundary of the domain, 
which is induced as the highly-conducting crust drags 
the frozen-in magnetic field, E — — [(fi + u>) x r] x £J, 
where CI = 2ire z /P represents rotation and u) represents 
shearing of the star. The outer boundary condition (at 
r ~ 3i?Lc) allows free escape of Poynting flux. The code 
has very low unphysical numerical diffusion or dissipa- 
tion on resolved scales, while a high-order spectral filter- 
ing procedure applies some controlled dissipation near 
the grid scale, which becomes the resistive length scale 
at which reconnection can occur. This is desirable, be- 
cause the magnetosphere should be very close to ideal, 
except in current sheets where reconnection is expected. 
For our purposes it is important only that resistivity be 
confined to current sheets, and that these can become 
as thin as possible — it is irrelevant that it is introduced 
with a filter rather than an explicit term in the equations 
of motion. 

In our axisymmetric simulations the magnetic axis is 
parallel to both the rotation and twist angular veloc- 
ity vectors, f2 and u:. Each simulation begins with the 
stationary solution for a rotating star with no twisting, 
lj = (Fig. la). It is found by setting a dipole magnetic 
field into rotation and allowing it t o relax to a steady 
state, as described in iParfrev et al.l (|2012ft . Then u> is 
superimposed, so that either a northern polar cap, or a 



TABLE 1 
Simulation parameters 



Run 




Twist Profile 


a 


N r X Ng 


A 


20 


Polar Cap 


6 


384 x 255 


B 


10 


Polar Cap 


12 


768 x 507 


C 


10 


Ring 


9-12 


768 x 507 



ring of northern latitudes, rotates slightly faster than the 
rest of the star. The polar cap, extending down to the 
colatitude pc , is modeled by the twisting profile 



>jj c 



w 



1 + exp [k (0 - pc )] ' 



and the ring profile is 

^ring(0) = — 



WO 

1 + exp {k [|0 — & c tr\ — A]} ' 



(4) 



(5) 



where 9 c t r is the colatitude of the center of the ring and A 
is its angular half-width. We set k = 50 and luo — fi/200; 
the twisting rate is approximately constant within the 
twisted region, and drops quickly to zero elsewhere. 

The twisted polar cap is bigger than the footprint of 
the open flux bundle of the (untwisted) rotating dipole, 

P c > 0rot, where rot ~ {r*/Rhc) 1/2 - The ratio of 
magnetic fluxes emerging through the caps 9 < pc and 
9 < rot defines a characteristic dimensionless parameter, 



f P c = sin 2 9 
/rot sin 2 rot 



pc 



> 1. 



(G) 



A similar dimensionless parameter can be used to label 
the flux surfaces that bound the ring twisting region. A 
ring defined by the flux surfaces a\ and a-i extends from 
d\ = 9 ctr — A to 2 = c tr + A, where sin 2 9\ = a\ sin 2 9 rot 
and sin 2 02 = 0-2 sin 2 ro t ; the magnetic flux through the 
twisted ring is (02 — <2i)/ TO t- 

Three simulations A, B, and C are detailed in Ta- 
ble 1, where N r and Ng are the numbers of grid points in 
the radial and meridional directions. In all simulations, 
the magnetosphere is twisted up to a maximum angle 
of Wntmax = 8 rads over 263 rotation periods. Fig. 1 
illustrates the behavior of run B. 

For runs A and B, the evolution proceeds as follows. 
As the twisting begins, toroidal magnetic field builds 
up on the closed field lines connected to the sheared 
cap (on open field lines the twist is emitted to infin- 
ity rather than accumulated). As the twist pumping 

ip is slow compared with the Alfven crossing time, the 
magnetosphere first evolves through a sequence of quasi- 
cquilibrium states. The increased magnetic pressure in 
the closed zone causes these field lines to expand out- 
wards, pushing more flux through the light cylinder and 
hence increasing the spindown torque on the star. The 
newly-opened field lines lose their twist — it is emitted to 
infinity — and establish the usual spindown-wind config- 
uration where the field changes sign across the equato- 
rial plane. The magnetic field discontinuity is supported 
by the equatorial current sheet, terminating at a Y-point 
which separates closed and open flux. The magnetic field 
remains almost reflection-symmetric about the equator 
(modulo differences in sign between oppositely directed 
flux) , because the twisting timescale is much longer than 
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(a) u>ot — (b) uot — 2.19 rads 
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(c) uot — 2.28 rads 
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Fig. 1. — The first two expansion and reconnection events from run B. Color shows toroidal magnetic field times the radial coordinate, 
in units where this quantity is initially ~ 1 at the light cylinder. Black curves show twenty-two poloidal field lines, equally spaced in 
poloidal flux function between / = 0.01/ ma x and (1/3) /max ! where /max = 2n/i/r ir ; one additional field line, that initially closes at the light 
cylinder, is shown in red. Axes are labeled in units of J?lc = 40r*. (a) The initial state (untwisted rotating star); (b) the last outward 
"breathing" motion before the first reconnection; (c) after the first reconnection, the inner dark blue region is the twisted reservoir, the 
almost white area is the cavity; (d) detonation — the reservoir expands explosively. 



the wave-crossing timescale on the closed field lines, al- 
lowing waves to distribute the twist nearly symmetri- 
cally. 

The expansion rate increases while the twist ip grows 
with constant rate w, and eventually the field lines ex- 
pand through the light cylinder faster than they adjust 
to a new spindown steady state; this sets up a pat- 
tern of magnetospheric "breathing," where the amount of 
open flux oscillates, via reconnection in the current sheet, 



with increasing amplitude and a period of ~ 7 Rlc/ c 
(Fig. lb). 

After a significant fraction of the polar cap flux has 
been opened, at tp ~ 2, the recently opened field lines be- 
come unstable to catastrophic reconnection and a large 
fraction of the open flux reconnects, bringing the Y-point 
within the light cylinder. The released magnetic energy is 
expelled in a plasmoid-fragmented outflow in the equato- 
rial plane. The Y-point quickly moves back to the light 
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1 2 3 4 5 6 7 
Footpoint angular displacement, cjot (radians) 

Fig. 2. — Spindown torque on the star (integrated angular mo- 
mentum flux through the stellar surface), in units of the torque on 
the untwisted rotating star, versus the relative angular displace- 
ment of the twisted region, ojo*. Runs A and B (polar caps) are in 
the upper panel, run C (ring profile) is in the lower panel. 



cylinder, opening up a cavity of zero toroidal field be- 
tween the low-lying strongly twisted region and the open 
flux bundle (Fig. fc). We will term the evolution from 
Fig. la to Fig. lc a "gradual" expansion and reconnec- 
tion event. 

The accumulated twist ip — ujQt remains intact on those 
lines, lying closer to the star, which never opened; let 
us call these field lines the "twisted reservoir." As the 
shear motion of the polar cap continues, toroidal field is 
added to both the twisted reservoir and the cavity. Be- 
fore a new significant ip builds up in the cavity, the over- 
twisted reservoir becomes unstable and expands outward 
explosively, and all the overlying field lines are briefly 
pushed through the light cylinder (Fig. Id). This process 
is rem iniscent of magnetic detonation (|Cowlev fe Artunl 
1997j). There is a narrow spike in the spindown torque, 
larger than in the previous event because more flux is 
opened. The enhanced open flux immediately recon- 
nects and closes, leaving a reservoir of lesser volume with 
ip ~ ujot and a bigger cavity with ip » 0. This is an "ex- 
plosive" expansion and reconnection event. 

This cycle of a twisted reservoir exploding and recon- 
necting is repeated several times as the shearing con- 
tinues. In the first few events successively more flux is 
opened, leading to higher spikes in spindown rate, and a 
deeper trough in torque immediately following reconnec- 
tion, since the twisted reservoir comprises ever fewer field 
lines. The spikes in spindown power become increasingly 
narrow as they become more energetic, with the most 
dramatic events having a duration of only a few light 
crossing times of the light cylinder (Fig. [2]). 

The time between events grows, because the twisted 
reservoir is shrinking. Eventually, the additional twist 
needed to detonate the reservoir becomes greater than 
~ 2 rads, and the cavity fills with enough toroidal field 



to initiate a gradual reconnection event before the reser- 
voir explodes. After this point the picture becomes more 
complicated, as there can be a twist-free cavity on top 
of a twisted reservoir comprising both flux that did not 
open in the last event, and flux that has never opened. 
Different amounts of flux will then expand and open de- 
pending on the history of the system. 

The heights of the torque spikes scale approximately 
as a 2 , as can be seen in Fig. [2j This scaling, which is 
confirmed by further simulations with a from 2 to 12, is 
expected. The spindown luminosity L is related to the 
magnetic flux through the light cylinder /lc by 



1 (he 

c I P 



(7) 



wher e P is the rotation period (jGoldreich fc Julian! 
1969). The opening of the twisted, inflated flux bundle 
increases /lc by nearly a factor of a, /lc ~ a/rot, which 
boosts L by a factor of a 2 . When averaged over the en- 
tire run time of the simulations, the spindown power in 
runs A and B is enhanced by factors of 5.45 and 9.98 
respectively; thus long-term average spindown is linear 
in a. 

Run C evolves slightly differently, because only deep 
closed flux is twisted. There is no gradual opening — 
all events are explosive, and require a larger twist an- 
gle than in the polar cap simulations before the accu- 
mulated magnetic pressure is great enough to overcome 
the restraining tension of the overlying untwisted field 
lines. The torque peaks scale as the square of the total 
opened flux (equation[7]) , which includes both the twisted 
ring and the overlying magnetosphere, because during an 
event the expanding twisted field lines push the overlying 
untwisted flux through the light cylinder. The average 
spindown power is enhanced by a factor of 2.54. 

3. DISCUSSION 

We find that the twisted magnetosphere undergoes 
a succession of opening and reconnection events, ac- 
companied by large increases in the spindown power 
(Fig. [2|). There are two main classes of field expan- 
sion and reconnection events. (1) The gradual kind in- 
volves smooth expansion of the magnetosphere through 
the light cylinder; it causes spindown to increase on the 
twisting timescale lu~ 1 until sudden reconnection restores 
the normal regime. (2) Explosive events are caused by 
the sudden breakout of flux under high magnetic pres- 
sure, whose expansion had been slowed by the overlying 
magnetosphere; these may be as brief as a few light cross- 
ing times of the light cylinder, which is comparable to the 
spin period of the star. In both cases the reconnection, 
and hence energy dissipation, phase is very short (as can 
be seen in the near- vertical downstrokes in Fig. [2]). 

The energy stored by twisting a polar cap 9 pc of a 
star with magnetic moment n and radius r* is given by 
(Bcl oborodovil2009D . 



M y sm 6 e pc 

24r 3 



4 x 10 4 V 2 3 ^ 



sin 6 6» r 



erg, (8) 



where ^33 = /i/10 33 G cm 3 and r* ss 10 6 cm. During an 
explosive event, this energy is partially released in a mag- 
netic outflow and partially radiated in a powerful flare. 
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Our simulations suggest that such events are accompa- 
nied by a brief increase in spindown torque by a factor 
as large as a 2 (see equation [7]), which leads to an abrupt 
increase in spin period. Then the fractional change in 
period is approximately given by 

AP 9 At , / cAt \ ■ 

a 2 a 2 — — P , (9 

where Pq — P/to is the spindown rate for the un- 
twisted magnetosphere and At is the duration of the 
torque increase by a 2 . A typical 9 pc ~ 0.3-0.6, which 
is sufficient for observed giant flares (equation |5J| , corre- 
sponds to a ~ 3 x 10 3 (equation EJ). SGR 1900+14 has 
Pq ~ 10" 10 s s" 1 , so a short-duration spike in spindown, 
with At comparable to Plc/c, can give AP/P ~ 10~ 4 
as observed for the August 1998 flare. This suggests that 
huge anti-glitches may be explained without recourse to 
sudden changes in the stellar interior. For the Decem- 
ber 2004 flare in SGR 1806-20, which was two orders of 
magnitude more energetic, one could expect even larger 
AP, which is not obser ved — the measured upper limit 
for AP/P is 5 x 10~ 6 (IWoods et al.l f2007h . Note that 
AP oc a 2 while E tw oc a 3 . Variations in a, /i, and twist 
geometry may lead to significant variations in AP. 
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particular kink instability of the inflated magnetosphere. 
We plan to investigate this possibility with future simu- 
lations. 
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